Vortices in a rotating Bose-Einstein condensate under extreme elongation in a 

harmonic plus quartic trap 
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The behaviors of a rapid rotating Bose-Einstein condensate under extreme elongation in a 2D 
anisotropic harmonic plus quartic trap are investigated. Due to the quartic trap, the system remains 
stable at high rotating velocity, > lo± (to± is the radial harmonic trap frequency), and vortex 
lattices form even in the absence of the repulsive s-wave interaction (g). When g is present, the 
interplay between g and the quartic trap potential can lead to rich vortex lattice transition states 
as a function of Q, to which vortex lattices vanish eventually at some higher Q. 
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I. INTRODUCTION 

In recent years vortices in a rotating Bose-Einstein con- 
densate (BEC) have been extensively studied theoreti- 
cally and experimentally. Rotating condensates are typ- 
ically confined in a harmonic potential such that large 
vortex arrays are obtained when angular velocity £1 is 
smaller than but approaching to the radial trap oscilla- 
tor frequency oj±. For instances in Refs. a single 
vortex state or a state of several hundreds of vortices have 
been successfully created. Theoretically, the behaviors of 
the rapid rotating BEC systems are usually studied based 
on the lowest Landau level (LLL) approximation (8l-[l3j. 
When angular velocity is greater than the radial trap os- 
cillator frequency, £1 > uj± , the system becomes unstable 
due to the ineffectiveness of confinement. 

In the experiment of Bretin et. al [14], a quartic 
trapping potential has been successfully created. Sev- 
eral groups then consider the type of trap potential 
V(x, y) = \muj\r\ + \ur\, where m is the atom mass, u 
is the strength of the quartic potential, and r\ = x 2 + y 2 
[l5l-[l9j. With this quartic potential, due to the effect of 
strong confinement, one thus can study the rotating Bose 
condensates at higher angular velocities, £1 > uj±. It is 
understood that the centrifugal force, varying as fl 2 r±, 
can always be compensated by the trapping force varying 
as — (muo^r^ + ur\). 

Rotating BECs have also been studied under an 
anisotropic harmonic trap for a quasi-two-dimensional 
(2D) system [20|, Moreover, by manipulating the 

anisotropic trap frequencies associated with the rotating 
velocity, quasi-2D rotating BEC under the anisotropic 
harmonic t rap can be reduced to the limit of extreme 
elongation [22-24]. In Ref. [22| the authors have studied 
the roton-maxon excitation spectrum in a rapid rotating 
weakly interacting BEC under extreme elongation in a 
2D anisotropic harmonic trap. When the interaction or 
the angular velocity increases, the system could undergo 
a second-order quantum transition to a state of a periodic 
structure - rows of vortices [22| . Similar transitions have 
also been investigated in Refs. [23|, [24[. A common find- 
ing is that the larger the interaction or angular velocity 



is, the more the vortex lattice row number is. 

Inspired by the works mentioned above [3, 1221424) , it 
is of particular interest to study how the rapid rotating 
BEC behaves under the extreme elongation if a quartic 
trapping potential is added. With this goal, this paper 
considers the following effective 2D trap potential 
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where, with the strength u, the anharmonic quartic po- 
tential acts on the x-axis only. By adding the quartic 
potential along the x axis only, the effect of the anhar- 
monic trap for the extreme elongated system can be seen 
in a more transparent way (effective ID problem). In 
fact, the kind of effective trap (pp) can be realized in a re- 
alistic experiment. One can consider a quasi-2D system 
that is confined in the xy plane by a harmonic trap of 
relatively large uo z . Now, consider that an extra isotropic 
2D quartic trap, V4 ~ (x 2 + z 2 ) 2 , is added and suppose 
that the rotating angular frequency is purely along the 
z direction. Since the system is already confined in the 
xy plane, any quartic trap effect in the z direction will 
simply enhance the confinement of the system in the xy 
plane. This means effectively, V4 ~ x A . 

With the quartic trapping potential (pQ) which is in- 
troduced to stabilize the system against decay at higher 
rotation velocities (Q), one can then study the system 
even for Vt > uo x . It will be shown that in the limit of 
no interaction, ground state forms no vortex lattice at 
smaller ft and when Q becomes larger, one row of vortex 
lattice will form. Vortex lattice will melt eventually at 
some larger Q. When the interaction is turned on, melt- 
ing of the vortex lattice at some higher Q is also found. 
Nevertheless, rich transition states occur for the vortex 
lattices. The transition states, arising due to the instabil- 
ity at a critical interaction strength or a critical rotation 
velocity, were first studied in Ref. [22| in an elongated sys- 
tem without the quartic potential. One will see that a 
similar instability leading to the formation of a periodic 
structure of vortex rows occurs in the current context. 
Of most interest, a parameter labeled by 7V C which cor- 
responds to number of terms involved in the /c-sum of 
the wavefunction is intimately related to the number of 
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vortex rows [22| . 

The paper is organized as follows. In Sec. Ull we out- 
line the GP energy functional essential for a rapid ro- 
tating BEC system under extreme elongation in a 2D 
anisotropic harmonic plus quartic trap. Details about 
how the ground-state wave function is obtained and how 
the total energy is minimized are given. In Sec. IIII| we 
study the properties of the system without the interac- 
tion. This section is devoted to fully understand the ef- 
fect of the quartic potential. It is shown that the be- 
haviors of the system with a quartic potential are quite 
different compared to those of no quartic potential. In 
Sec. IIV( effect of the s-wave interaction g is studied, 
which is treated by the perturbation method. It will 
be shown that the interplay between g and the quartic 
trap potential can lead to rich vortex lattice transition 
states as a function of O, to which vortex lattices vanish 
eventually at some higher Q. Sec. [V] is a brief conclusion. 



II. GP ENERGY FUNCTIONAL 

The GP energy functional of a rapid rotating, quasi-2D 
BEC system can be given by 



y direction. This is the case of most interest in the cur- 
rent context (effective ID case). Nevertheless, for prac- 
tical reason, we will assume that the system is confined 
in a large rectangular box of side length L. When u = 
(no quartic trap), as mentioned before, the system will 
become unstable for Q > ujo unless the s-wave interac- 
tion is attractive (#2D < 0). However when u ^ 0, the 
system remains stabilized for Q > uo and #2D > 0. 

To solve the wavefunction of the system, we shall use 
an approach similar to that used in Refs. [22, 23|. Firstly, 
we apply the gauge transformation, = \&e~ l7nnxy / h , 
and consequently the energy functional (j4]) becomes 



E 



J dxdy^l - 
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where the trap potential V(x,y) was given in Eq. (pp), 
= ^o( x ,y) is the wave function, #2D = 
4nh 2 a s /(ma z \ /r 27r) is the effective 2D s-wave interaction 
with a s the s-wave scattering length and a z = y/h/muj z 
associated with the harmonic trap frequency u z [25j, CI 
is the rotation angular velocity, and L is the angular mo- 
mentum. When ft = Qz, the rotating term becomes 



Eq. (j2j) can also be written as 



d_ 

dy 



d_ 

1 dx 



^0- 



E = / dxdy% 



(—ihV + mr± x fi) 
2m 



Ves(x,y) 



9 ~f\*o\ 2 



(3) 



(4) 



where = (x,y) and V G s(x,y) = (m/2)[{uj 2 - Q 2 )x 2 + 
(uo 2 — Q 2 )y 2 ] + (u/4)x A . For the part of the anisotropic 
harmonic trap, trap frequencies can be written as uj 2 = 



ujq + A 2 and uo 2 = uu Q 
effective potential becomes 



uo 2 — A 2 . If one sets A 2 = uo 2 — Q 2 , the 



V e ff(x, y) = m(ujQ - n 2 )x 2 + -i 



(5) 



leaving no y dependence in it. This means that the sys- 
tem will reach a state under extreme elongation along the 
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with Veff(x,y) = m(uj 2 + £l 2 )x 2 + (u/4)x 4 . Secondly 
the s-wave interaction #2D is assumed to be weak such 
that wavefunctions corresponding to no interaction case 
(#2D = 0) can be used as the basis functions. Thirdly, 
the effect of the s-wave interaction will be reinstalled and 
treated by the perturbation method. Thus it is proposed 
that the wave function ^ is a linear combination of the 
lowest level eigenfunction(s), associated with energy 
functional (j6]) with # 2 d = 0. That is, 



9 = VNj2Ck<f> k , 



where N is the total number of particles and it requires 
that \Ck\ 2 = 1 for completeness. Here k = 2tt£/L 
with i an integer. Since L is large, k can be treated in the 
continuous limit. When the quartic trap vanishes (u = 
0), <pk corresponds to LLL. Furthermore, eigenfunctions 
4>k are assumed to take the separate form 



(j) k = Ae iky X k (x) , 



(8) 



where A = 1/y/L is the normalization constant. Sub- 
stituting Eqs. (|7J) and © into Eq. (j6j) with g = 0, one 
obtains the one-dimensional (dimensionless) differential 
equation for Xk (x), 



dx 2 



2Q 2 



-kx + Xx 4 + k 1 



: £kXk, (9) 



where is the eigenenergy. Here A = uh/2m 2 ujQ and 
all the le ngths a nd energies are scaled by the unit length 
ao = yh/muoo and unit energy Eq = hwo/2. When 
A = 0, one has analytic solutions for Eq. (j9j) [23]. But 
for A ^ 0, no analytic solution is available and numer- 
ical computation is needed. Here we shall apply the fi- 
nite difference method to find solutions of the normalized 
function Xk- 
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FIG. 1. (Color online) In frame (a)-(c), lowest eigenenergy 
of Eq. J9}, is plotted as the function of k for three cases: 
LJo/Q — 1.0,0.99 and 0.98 respectively. Frame (d) plots N c 
(see text) as the function of ujo/ft. In each frame, the four 
lines correspond to A = 0.0, 0.005, 0.01 , and 0. 03 respectively. 
For the units, Eq — huo/2 and ao — ^/h/mujo. 



With Eqs. 0-©, the full GP energy functional @ 
becomes 



E ng 
2 

k ki ,/C2 ,/C3 ,/C4 

X / dxXk 1 Xk 2 Xk 3 Xk A $k 3 +k A ,k 1 +k 2 , 



(10) 



where n = N/(L/ao), g = 4v27ra s / 'a z represents the new 
scale of #2D, and N c is the number of terms used in the 
linear combination. 7V C is determined upon the condi- 
tion that energy functional is minimized. The first part 
of the GP functional ([TO]) involves the kinetic, harmonic, 
anharmonic (quartic), and rotational energies. While the 
second part corresponds to the interaction energy ob- 
tained by the perturbation manner. It is worth noting 
in Eq. ((9]) that the function Xk(%) — X-k(-x) and the 
energy Sk = S-k- Due to the symmetry of Sk = it is 
useful to write Eq. ([7]) as 



. (11) 



With Eq. (fTT]h it's evident that there are two cases for 
N c . One is N c = 2j + 1 if C ^ and another is N c = 2j 
if Co = 0. 



III. NONINTERACTING SYSTEM 

In this section, we consider the ground-state properties 
of a noninteracting system (g = 0). It will be shown for 



this noninteracting system that when A / 0, the prop- 
erties of the system are quite different from those of the 
A = case. It is useful to first review the results for 
the A = case [23|. When A = 0, one has the analytic 
solution for energy Sk in Eq. (|9]) 
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The results of Eq. (fl2j) can be divided into three folds, 
(i) When ft < cjo, minimum e k occurs at k = 0. Con- 
sequently the ground-state wavefunction of the system is 
just = VNCo^o (N c = 1). (ii) When ft = u , e° k is 
the same for all fc's. In this highly degenerate case, the 
ground-state wavefunction should include all eigenstates 
(N c = oo). (hi) When ft > cjq, minimum e k occurs at 
k — >• oo. It indicates that the system is unstable. An im- 
portant consequence of the g = A = case is that vortex 
(lattice) will not form regardless of the ratio of uo/ft. 

When A 7^ 0, in contrast, behaviors of the system can 
be very different. We first consider the case when A is 
small such that wavefunction of the A = case can be 
used to obtain the A correction perturbatively. The lead- 
ing correction of the energy is found to be (ek — s k -\-Ae k ) 



Ae° k = X(A 



A 2 k 2 +A A k A ), 



(13) 



where A = 3/(8iJ), A 2 = 3ft 2 /(V2uj 2 R 5 ^ 2 ), and A 4 = 
ft 4 /(u%R 4 ) with R = 1 + ft 2 /u 2 . Note that all three A { 
are positive, so the correction Ae k > 0. Moreover, it is 
found that the behaviors of the wavefunction \I/ can be 
divided into two branches. When ft < (1 + 3A/8)cjo, min- 
imum Sk occurs at k = 0. Consequently \I> = y/NCofo 
(N c = 1). In contrast when ft > (1 + 3A/8)o;o, ^ = 
V^(Cfci0fci+C'-fci0-fci) (N c = 2), where k± corresponds 
to where minimum Sk is. Comparing these small but fi- 
nite A results to those of A = discussed before, one sees 
clearly that the general behaviors are quite different be- 
tween the A = and A ^ cases, especially in the regime 

ft > UJQ. 

When A is not too small such that perturbation ap- 
proach is no longer valid, one needs to solve Eq. (j9j) nu- 
merically for Sk- In Fig. QJa)-(c) with uo/ft = 1.0, 0.99, 
and 0.98 respectively, we numerically solve and plot Sk 
as the function of k. Three finite A cases (A = 0.005, 
0.01, and 0.03) are considered. The A = case is in- 
cluded for comparison. Based on the experimental data 
given in Ref. [14] that w x ~ 2tt x 64.8 Hz corresponds to 
A ~ 0.001, it is estimated that w x ~ 2ir x 38.5, 2tt x 30, 
and 2tt x 21 Hz will correspond to A = 0.005, 0.01, 
and 0.03 cases respectively. As shown in Fig. [TJa) for 
ujo/ft = 1.0, k\ (corresponding to minimum Sk) is zero for 
all three finite A cases. This is in big contrast to the A = 
case where e k is k independent [see Eq. ([T2]) ]. In the case 
of ujo/ft = 0.99, k\ is finite for A = 0.01 and 0.005 cases 
but is zero for the A = 0.03 case [see Fig. QJb)]. While 
the corresponding e k becomes lesser when k is larger and 
no minimum-energy state is found. Finally for the case 



(a)coo/Q=1.0 



(b) gVQ=0.995 




FIG. 2. (Color online) Normalized atom number density dis- 
tribution, |^| 2 /n, is plotted in the xy plane for (a) uo/Cl = 
1.0, (b) 0.995, (c) 0.993, and (d) 0.99. Here A is 0.005 for all 
frames. As shown, single vortex line exhibits in frame (b) and 
(c) only, x and y axes are in units of clq = ^h/muoo. 



of ujq/CI = 0.98 [see Fig. ffl^c)], k\ is finite for three A ^ 
cases. Similar to the uo/Cl = 0.99 case, will become 
lesser when k is larger and no minimum-energy state is 
found. Thus for the larger A (= 0.03) case, vortex forms 
only when the angular velocity (Cl) is larger. Moreover 
with the same A, the value of k\ (corresponding to mini- 
mum £fc) is larger when Cl is larger. This means that the 
number of vortices per unit length (in the elongated di- 
rection) will be larger when Cl is larger. In Fig. QJd), N c 
is plotted as the function of ujo/CI for different A cases. 

In practice, it is useful to compare the energy scale of 
the (dimensionless) quartic coupling A to the harmonic 
trap frequency uo x or uo y . When Cl ~ cjo, = 2 — 

(uo/Cl) 2 ~ 1. Since the results discussed above are for 
A < 0.03, the effect of the quartic trap is seen to be quite 
drastic. 

In Fig. [2j normalized atom number density distribu- 
tions, |^| 2 /n, are plotted in the xy plane. Four cases are 
considered, namely uoq/CI = 1.0, 0.995, 0.993, and 0.99 
with A = 0.005 in all cases. In the case of uoq/CI = 1.0 
[Fig. [2fa)], no vortex is formed. When the rotation ve- 
locity Cl is increased, single vortex line is seen to occur 
for both the uo^/Cl = 0.995 and 0.993 cases [Fig. E{b) 
and (c)]. It is noted, however, that number of vortices 
per unit length is larger for the ujo/CI = 0.993 case as 
compared to that of the uoq/CI = 0.995 case. When Cl is 
further increased to above a critical value, uo/Cl c < 0.99, 
vortex state disappears (melts) due to the large centrifu- 
gal force [see Fig[2jd)]. At this large Cl case, atoms are 
pushed to two sides along the elongated potential well 
and no atoms are left in the center. 

The most important results obtained in this section 
for the noninteracting system (g = 0) are summarized as 
the following, (i) The quartic trap can make the system 
remain stable at higher rotation velocity (Cl > ojq). (ii) 
Single vortex line can exhibit at some Cl regimes, (iii) The 
single-line vortex lattice will vanish eventually at some 
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FIG. 3. N c is plotted as the function of uo/Cl for six cases: (a) 
(ri0, A) = (1,0.03), (b) (1,0.01), (c) (1,0.005), (d) (0.5,0.03), 
(e) (0.5,0.01) and (f) (0.5,0.005). Transition states occur in 
all frames except (d). 



higher Cl to which atoms are push to two sides along the 
elongated potential well. In next section, the effect of 
interaction (g ^ 0) will be discussed in details. 



IV. EFFECT OF INTERACTION 

In this section, the effect of interaction (g) is studied 
by the perturbation approach. Similar to previous sec- 
tion, we consider also the three finite A cases: A = 0.03, 
0.01, and 0.005. When the energy associated with the 
quartic trap is smaller than the energy associated with 
the harmonic trap and when Cl ~ co>o, the energy differ- 
ence between lowest and first excited band of Sk is about 
2/kJo- In our dimensionless scale, 2hujo = 4. Thus the 
perturbative approach is valid as long as 



||^| 2 <4 or 



1*1 



ng 



«8. 



(14) 



The value of |*| 2 /n is roughly 0.4 (see Figs. [2]and[5j). 
Thus when ng <C 20, the perturbation approach is valid. 
Hence the values of ng = 0.5 and 1.0 considered here 
should be safely valid for the perturbative approach. 
Typically n is more than 1 in a 2D system [22|. If 
N = 10 4 and L/a Q = 100, the value of n is about 100. 
For ng = 1, it corresponds to a s /a z ~ 10 -3 . In Fig. [3j 
the number of N c [terms involved in the summation ([TO]) ] 
is determined and plotted as the function of uo/Cl for 
six cases, namely (ng, A) = (1,0.03), (1,0.01), (1,0.005), 
(0.5,0.03), (0.5,0.01), and (0.5,0.005). It turns out that 
the results of A^ c are quite fascinating in terms of the 
change of Cl and the interplay between the strength of 
the interaction and the strength of the quartic potential. 
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(aWQ=1.02,Nc=2 



(b)<Bo/n=1.0,Nc=3 
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FIG. 4. N c is plotted as the function of ng for three cases: (a) 
ujo/tt = 0.9815, (b) 0.99, and (c) 1.0. A = 0.005 in all three 
frames. 



One common feature of all frames in Fig. [3] is that N c 
shifts from 1 (at lower Q) to 2 eventually (at higher Q). 
In between of the 1 to 2 period, transition states occur 
for a large Q span. For example, N c > 3 transition states 
can occur in all frames except in frame (d). Besides, 
N c increases in a trend as Q increases for the transition 
states. The reason why N c all shift from 1 (at lower fi) 
to 2 eventually (at higher Q) can be understood as fol- 
lows. One recalls the N c results for g = in Fig. GJd), 
where all three finite A cases are observed to exhibit the 
N c = 1 — >> 2 transition (at different critical Q though). 

When the interaction g is turned on and still valid in 
the perturbative regime, g will only play little role in the 
large £1 limit. It means that when £1 is large enough, 
transition states, which arise due to the effect of g, will 
disappear. As an extreme example shown in Fig. Efd), 
because g is relatively small (ng = 0.5) and the quartic 
coupling is relatively large (A = 0.03), consequently no 
transition state occurs. When A is reduced [see Fig. [3fe) 
and (f)], or when g is increased [see Fig. [3ja)], transition 
states will occur. 

Moreover, it is found that the critical value of Q to 
which N c changes from 1 to 2 at the lower ft side is 
smaller when A is smaller (if ng is fixed) or when ng 
is larger (if A is fixed). It implies that the system will 
enter the vortex state earlier if g is relatively larger or 
A is relatively smaller. Furthermore, with the same A, 
transition states will sustain for a larger span of Q if g 
is larger and can go up to a higher value of N c . The 
latter means that vortex lattice can have a higher row 
number. One sees in Fig. [3](c) that N c can go up to 5, 
although for a small period of Q. Fig.HJplots number N c 
as the function of ng with A = 0.005 and ujo/Q = 0.9815, 
0.99, and 1.0 respectively. Basically Figs. [3] and |4] show 
how the different regimes behave as the change of the 
interaction, quartic trap strength, and rotation velocity. 

The normalized density profile in the xy plane, |\I/| 2 /n, 
is shown in Fig. [5]for (ng,\) = (1,0.005) [corresponding 
to the case in Fig.[3fc)]. As seen in Fig. [3jc), the change 
of N c in the transition states is quite rich for this case. Six 
angular velocities are studied, namely uo/Q = 1.02, 1.0, 
0.984, 0.9815, 0.98 and 0.9785 respectively for Fig.G^a)- 
(f). With these values of ^, N c corresponds to 2, 3, 4, 




FIG. 5. (Color online) Normalized atom number density dis- 
tribution, |^| 2 /n, is plotted in the xy plane for (a) uo/Q = 
1.02, (b) 1.0, (c) 0.984, (d) 0.9815, (e) 0.98, and (f) 0.9785. 
Here ng = 1 and A = 0.005 for all frames [corresponding to 
the case in Fig. C2c)]. Vortex lattices with the transitions of 
row number and lattice constant ar e obser ved as Q changes. 
x and y axes are in units of ao = ^h/mujQ. 



5, 4, and 2 respectively. One sees in Fig. E^a)-(d) that 
when Q is increased from ujo/Q = 1.02 to 0.9815, atoms 
are pushed to the two sides and the number of vortices 
becomes more and more. The row number of vortex line 
also increases from 1 to 4. While in Fig. Efe), the vor- 
tex row number is reduced to 3 (N c = 4) again. In the 
case of Fig. 0(f), although N c is reduced to 2, but vortex 
lattice vanishes (melts) due to the large centrifugal force. 
Similar vortex lattice melting transition (atoms are com- 
pletely pushed to the two sides) at large ft has already 
been seen in the previous section of no interaction. 



V. CONCLUSIONS 

This paper investigates the effect of a quartic potential 
on a fast rotating BEC system under the extreme elonga- 
tion. In contrast to the harmonic trap alone case where 
system is unstable when the angular velocity Q is larger 
than the radial trap oscillator frequency c<Jo, the quartic 
trap can lead the system to remain stable at higher rota- 
tion velocity (Q > ljo). The interplay between the weak 
s-wave interaction and the quartic trap potential can re- 
sult rich vortex lattice transition states as a function of 



6 



Q. At large f2, atoms are eventually push to the two sides also acknowledge the support from NCTS, Taiwan, 
along the elongated potential well. 
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